#include "Path_Hop.h"
Path_Hop::~Path_Hop()
{
    delete p_interval_t, p_interval_nt, p_intl;
    free(dfs_array);
    for (int i = 0; i < node_count; i++)
    {
		free(tree_interval[i]);
		free(utc[i]);
    }
    free(tree_interval);
    free(utc);
}

int Path_Hop::single_edge_connection()
{
    int orphan_cnt = 0;
    map<int, set<int> >::iterator it_m1, it_m2, it_m;
    set<int>::iterator it_s, it_s1;
    int a, b;
    bool single_edge;
    for (it_m = descendants_list.begin(); it_m != descendants_list.end(); it_m++)
    {
		single_edge = true;
		a = it_m->first;
		for (int i = 0; i < node_count; i++)
		{
			if (!utc[i][a])
			{
				single_edge = false;
				break;
			}
		}
		if (!single_edge) continue;

		for (it_s = it_m->second.begin(); it_s != it_m->second.end(); it_s++)
		{
			b = *it_s;
			if (utc[a][b]) continue; // already covered
			single_edge = true;
			for (int i = 0; i < node_count; i++)
			{
				if (!utc[b][i])
				{
					single_edge = false;
					break;
				}
			}
			if (!single_edge) continue;

			for (int c = 0; c < node_count; c++)
			{
				//if (c == a ) continue;
				if (!utc[a][c] && !utc[c][b])
				{
					single_edge = false;
					break;
				}
			}
			if (single_edge)
				orphan_cnt++;
		}
    }
    /*
      for (int a = 0; a < node_count; a++)
      {
      for(int b = 0; b < node_count; b++)
      {
      if (!utc[a][b])
      {
      single_edge = true;
      for (int c = 0; c < node_count; c++)
      {
      if ( c == a) continue;
      if (!utc[a][c] && !utc[c][b])
      {
      single_edge = false;
      break;
      }
      }
      if (single_edge)
      orphan_cnt ++;
      }
      }
      }
    */
    /*
      for (it_m = ancestors_list.begin(); it_m != ancestors_list.end(); it_m++)
      {
      for (it_s = it_m->second.begin(); it_s != it_m ->second.end(); it_s++)
      {
      if (ancestors_list.find(*it_s) == ancestors_list.end())
      utc[*it_s][it_m->first] = false;
      }
      }
    */
    //cout << "#Single Edge Connections:" << orphan_cnt << endl;
    return orphan_cnt;
}

/*
  reduce the graph into a minimal equivalent graph
*/
void Path_Hop::reduce_graph(Graph& g)
{
    graph_traits<Graph>::in_edge_iterator ei, eiend;
    graph_traits<Graph>::in_edge_iterator ei2, ei2end;
    set<size_t> parents;
    set<size_t>::iterator p_it;
    vector<Vertex>::reverse_iterator ri;
    Vertex v, u;
    Vertex u_tc, ancestor;
    Edge e;
    bool has_edge;
    size_t pid;

    for (ri = topo.rbegin(); ri != topo.rend(); ++ri)
    {
		v = *ri;
		g[v].x = 0; // used for finding spanning tree
		parents.clear();
		/* check redundant edges */
		// get all parents
		for (tie(ei, eiend) = in_edges(v, g); ei != eiend; ++ei)
		{
			u = source(*ei, g); // get each parent
			parents.insert(g[u].id); // store the id of parent
		}
		// process each parent
		for (p_it = parents.begin(); p_it != parents.end(); ++p_it)
		{
			pid = *p_it;
			tie(e, has_edge) = edge(vertex(pid, g), v, g);
			if (!has_edge)
			{
				continue;
			}
			u_tc = vertex(pid, tc);
			// get u's ancestor
			for (tie(ei2, ei2end) = in_edges(u_tc, tc); ei2 != ei2end; ++ei2)
			{
				ancestor = source(*ei2, tc); // get each ancestor
				assert(ancestor != u_tc);

				// remove redundant edge
				if (parents.find(tc[ancestor].id) != parents.end())
				{
					remove_edge(vertex(tc[ancestor].id, g), v, g);
				}
			}
		}
    }
}

/////////////////////////////////////////////////////////////////////////////////

void Path_Hop::topo_sort()
{
    //compute the transitive closure
    std::vector<Vertex> to_tc_vec(num_vertices(gra));
    VertexIndexMap index_map = get(&VertexProp::id, gra);
    iterator_property_map < Vertex *, VertexIndexMap, Vertex, Vertex&> g_to_tc_map(&to_tc_vec[0], index_map);
    transitive_closure(gra, tc, g_to_tc_map, index_map);
    graph_traits<Graph>::vertex_iterator vi, vi_end, v2, v2e;
    for (tie(vi, vi_end) = vertices(gra), tie(v2, v2e) = vertices(tc); vi != vi_end; ++vi, ++v2)
    {
		tc[*v2].id = gra[*vi].id; // copy ids
    }
    // topological sort
    try
    {
		topological_sort(gra, std::back_inserter(topo), vertex_index_map(get(&VertexProp::id, gra)));
    }
    catch (not_a_dag)
    {
		cout << "Error: Not a DAG for topological sort" << endl;
		assert(false);
    }
}

/*
  check if any nodes in a node's Lin has A-D relationship (redundancy), if so, remove A. keep D.
*/
void Path_Hop::remove_redundant_Lin_Lout_in_vertex()
{
    int Lin_cnt = 0, Lout_cnt = 0;
    map<int, set<int> >::iterator it_m;
    set<int>::iterator it_s, it_s1;
    set<int> overlap_set;
    //check Lout is not necessary
    for (it_m = Lout.begin(); it_m != Lout.end(); it_m++)
    {
		for (it_s = it_m->second.begin(); it_s != it_m->second.end(); it_s++)
		{
			overlap_set.clear();
			for (it_s1 = it_m->second.begin(); it_s1 != it_m->second.end(); it_s1++)
			{
				if (*it_s != *it_s1 && is_descendant(*it_s, *it_s1))
				{
					overlap_set.insert(*it_s1);
				}
			}
			for (it_s1 = overlap_set.begin(); it_s1 != overlap_set.end(); it_s1++)
			{
				it_m->second.erase(*it_s1);
				Lout_cnt++;
			}
		}
    }
    //cout << "Redundant Lout in each single node: " << Lout_cnt << endl;
    for (it_m = Lin.begin(); it_m != Lin.end(); it_m++)
    {
		for (it_s = it_m->second.begin(); it_s != it_m->second.end(); it_s++)
		{
			overlap_set.clear();
			for (it_s1 = it_m->second.begin(); it_s1 != it_m->second.end(); it_s1++)
			{
				if (*it_s != *it_s1 && is_descendant(*it_s1, *it_s))
				{
					overlap_set.insert(*it_s1);
				}
			}
			for (it_s1 = overlap_set.begin(); it_s1 != overlap_set.end(); it_s1++)
			{
				it_m->second.erase(*it_s1);
				Lin_cnt++;
			}
		}
    }
    //cout << "Redundant Lin in each single node: " << Lin_cnt << endl;
}

void Path_Hop::remove_redundant_Lin_Lout_in_tree()
{
    //Remove the redundancy in Lout, we start from the most bottom leaf nodes and traverse upwards.
    //We maintain a set of all Lout when traversing upwards and remove any redundant Lout if the Lout of a node is already contained in that set.
    int lout_cnt_r = 0, lin_cnt_r = 0;
    map<int, set<int> >::iterator it_in, it_out, it_m;
    set<int>::iterator p;
    vector<int>::iterator it;
    set<int> lout_set, lin_set;
    int parent;
    /*
      for (int i = 0; i < node_count; i++)
      {
      it_m = Lout.find(i);
      if (it_m == Lout.end()) continue;
      lout_set.clear();
      int begin = tree_interval[i][0];
      int interval = tree_interval[i][1] - tree_interval[i][0];
      for (int i = begin + 1; i <= begin + interval; i++)
      {
      it_out = Lout.find(dfs_array[i]);
      if (it_out != Lout.end())
      {
      lout_set.insert(it_out->second.begin(), it_out->second.end());
      }
      }
      for (p = it_m->second.begin(); p != it_m->second.end(); p++)
      {
      if (lout_set.find(*p) != lout_set.end())
      {
      it_m->second.erase(p);
      lout_cnt_r++;
      }
      }
      }
      cout << "Redundant Lin in whole tree: " << lout_cnt_r << endl;
    */

    for (int i = 0; i < node_count; i++)
    {
		it_m = Lin.find(i);
		if (it_m == Lin.end()) continue;
		lin_set.clear();
		parent = tree_interval[i][2];
		while (parent != -1)
		{
			it_in = Lin.find(parent);
			if (it_in != Lin.end())
			{
				lin_set.insert(it_in->second.begin(), it_in->second.end());
			}
			parent = tree_interval[parent][2];
		}
		for (p = it_m->second.begin(); p != it_m->second.end(); p++)
		{
			if (lin_set.find(*p) != lin_set.end())
			{
				it_m->second.erase(p);
				lin_cnt_r++;
			}
		}
    }
    //cout << "Redundant Lin in whole tree: " << lin_cnt_r << endl;
}

// for correctness check

void Path_Hop::post_process()
{
    p_intl = new Interval(gra);
    p_intl->createLabels(false);
    /*
      int tc = 0;
      for (int i = 0; i < node_count; i++)
      {
      for (int j = 0; j < node_count; j++)
      {
      if (reach(i, j)) tc++;
      }
      }
      cout << "#TC Size is :" << tc << endl;
    */
}

//convert the Lin Lout form set to vector, to speed up

void Path_Hop::set_2_vector()
{
    Lin_vec = vector<vector<int> >(node_count, vector<int>());
    Lout_vec = vector<vector<int> >(node_count, vector<int>());
    map<int, set<int> >::iterator it;
    for (int i = 0; i < node_count; i++)
    {
		it = Lin.find(i);
		if (it != Lin.end())
		{
			Lin_vec[i].insert(Lin_vec[i].begin(), it->second.begin(), it->second.end());
		}
		it = Lout.find(i);
		if (it != Lout.end())
		{
			Lout_vec[i].insert(Lout_vec[i].begin(), it->second.begin(), it->second.end());
		}
    }
}

int Path_Hop::num_special_nodes()
{
    int cnt = 0;
    vector<Vertex>::iterator it;
    for (it = topo.begin(); it != topo.end(); it++)
    {
		if (has_non_tree_edge(*it))
			cnt++;
    }
    cout << "special nodes: " << cnt << " , ratio: " << (float) cnt / node_count << endl;
}

/************************
/*test if a vetex has non tree edges or not, if not, then it is an ordinary node,i.e. tree nodes
/* result = 0 : has neither incomming nor out-going non-tree edges
/* result = 1 : has only out_going non-tree edges;
/* result = 2 : has only in_coming non-tree edges;
/* result = 3 : has both in_comming and out-going non-tree edges 
*************************/
int Path_Hop::has_non_tree_edge(Vertex v)
{
    int result = 0;
    graph_traits<Graph>::out_edge_iterator ei, eiend;
    for (tie(ei, eiend) = out_edges(v, gra); ei != eiend; ++ei)
    {
		if (!gra[*ei].tree)
		{
			result += 1;
			break;
		}
    }
    if (in_degree(v, gra) > 1)
		result += 2;
    return result;
}

int Path_Hop::has_non_tree_edge(int nid)
{
    Vertex v = vertex(nid, gra);
    return has_non_tree_edge(v);
}

void Path_Hop::init_tree_interval()
{
    tree_interval = (int**) malloc(node_count * sizeof (int*));
    set<IntPair>::iterator it;
    Vertex v, p;
    graph_traits<Graph>::in_edge_iterator ine, ine_end;
    // tree interval : start, end, parent_id
    for (int i = 0; i < node_count; i++)
    {
		tree_interval[i] = (int*) malloc(3 * sizeof (int));
		it = p_interval_t->intervals.find(i)->second.begin();
		tree_interval[i][0] = it->first;
		tree_interval[i][1] = it->second;
		tree_interval[i][2] = -1; // -1 is for root which has no parent
		v = vertex(i, gra);
		for (tie(ine, ine_end) = in_edges(v, gra); ine != ine_end; ine++)
		{
			if (gra[*ine].tree)
			{
				p = source(*ine, gra);
				tree_interval[i][2] = gra[p].id;
				break;
			}
		}
    }
}

/*
  jump parent pointer
*/
void Path_Hop::init_parent_pointer()
{
    int parent;
    for (int i = 0; i < node_count; i++)
    {
		parent = tree_interval[i][2];
		while (parent != -1)
		{
			if (Lin.find(parent) != Lin.end())
			{
				tree_interval[i][2] = parent;
				break;
			}
			parent = tree_interval[parent][2];
		}
		tree_interval[i][2] = parent;
		//cout << "Node " << i << " [" << tree_interval[i][0] << ", " << tree_interval[i][1] << "] " << tree_interval[i][2] << endl;
    }
}

/*
  check if v is u's descendant; i.e. u --> v;
*/
bool Path_Hop::is_descendant(Vertex u, Vertex v)
{
    return ( gra[u].begin < gra[v].begin && gra[u].end >= gra[v].end);
}

/*
  Check if v is u's descendant; i.e. u --> v; 
  Overloaded function;
*/
bool Path_Hop::is_descendant(size_t top_nid, size_t bottom_nid)
{
    return tree_interval[top_nid][0] < tree_interval[bottom_nid][0] && tree_interval[top_nid][1] >= tree_interval[bottom_nid][1];
}

/*
  Test if u's non-tree intervals has v's tree interval  (containment, not exactly equals)
*/
bool Path_Hop::has_interval(Vertex u, Vertex v)
{
    set<IntPair, pair_comp>::iterator it;
    map<int, set<IntPair, pair_comp> >::iterator iter;
    iter = p_interval_nt->intervals.find(gra[u].id);
    int begin = gra[v].begin;
    int end = gra[v].end;
    for (it = iter->second.begin(); it != iter->second.end(); it++)
    {
		if (it->first == begin && it->second == end)
		{
			return true;
		}
    }
    return false;
}


int Path_Hop::interval_2_nid(int begin, int end)
{
    map<int, PairSet >::iterator it_t;
    for (it_t = p_interval_t->intervals.begin(); it_t != p_interval_t->intervals.end(); it_t++)
    {
		if (it_t->second.begin()->first == begin && it_t->second.begin()->second == end)
		{
			return it_t->first;
		}
    }
    return -1;
}

int Path_Hop::begin_2_nid(int begin)
{
    for (int i = 0; i < node_count; i++)
    {
		if (tree_interval[i][0] == begin)
			return i;
    }
}

void Path_Hop::construct_descendants_sort_list()
{
    std::vector<Vertex>::reverse_iterator it, it_d, it_a;
    set<nid_vol_pair, vol_Compare> sorted_descendants;
    for (it = topo.rbegin(); it != topo.rend(); ++it)
    {
		//we only consider "non-tree nodes" here
		if (has_non_tree_edge(*it) >= 2) // it should has in-vol
		{
			//construct descendants sorted list
			sorted_descendants.clear();
			for (it_d = it + 1; it_d != topo.rend(); ++it_d)
			{
				if (((has_non_tree_edge(*it_d) % 2) == 1) && is_descendant(*it, *it_d)) // it should has out-vol as well as it is an tree descentdant
					//if( is_descendant(*it, *it_d) )
				{
					int out_vol = p_interval_nt->intervals.find(gra[*it_d].id)->second.size();
					if (out_vol > 0)
					{
						sorted_descendants.insert(make_pair(gra[*it_d].id, out_vol));
					}
				}
			}
			if (sorted_descendants.size() > 0)
				descendants_sort_list.insert(make_pair(gra[*it].id, sorted_descendants));
		}
    }
}

int Path_Hop::get_first_of_desc_sort_list(int v)
{
    map<int, set<nid_vol_pair, vol_Compare> >::iterator it_d;
    it_d = descendants_sort_list.find(v);
    if (it_d == descendants_sort_list.end())
		return -1;
    int bottom_v = it_d->second.begin()->first;
    int fake_out_vol = it_d->second.begin()->second;
    int out_vol = calc_out_vol(bottom_v); // real time value
    while (fake_out_vol != out_vol) //   fake > real
    {
		//remove old value
		it_d->second.erase(it_d->second.begin());
		//insert new value
		if (out_vol > 0)
			it_d->second.insert(make_pair(bottom_v, out_vol));
		if (it_d->second.size() > 0)
		{
			bottom_v = it_d->second.begin()->first;
			fake_out_vol = it_d->second.begin()->second;
			out_vol = calc_out_vol(bottom_v);
		}
		else
		{
			descendants_sort_list.erase(it_d);
			return -1; // -1 means not exists
		}
    }
    return bottom_v;
}

void Path_Hop::print_path_hop_labels(int node_count)
{
    cout << "print path-hop labels..." << endl;
    map<int, set<int> >::iterator it_in, it_out;
    set<int>::iterator it;
    int i = 0;
    while (i < node_count)
    {
		it_in = Lin.find(i);
		it_out = Lout.find(i);
		cout << i << ": ";
		if (it_in != Lin.end())
		{
			cout << "Lin(";
			for (it = it_in->second.begin(); it != it_in->second.end(); it++)
			{
				cout << *it << ", ";
			}
			cout << "); ";
		}
		if (it_out != Lout.end())
		{
			cout << "Lout(";
			for (it = it_out->second.begin(); it != it_out->second.end(); it++)
			{
				cout << *it << ", ";
			}
			cout << ")";
		}
		i++;
		cout << endl;
    }
}

int Path_Hop::path_hop_index_size()
{
    map<int, set<int> >::iterator it;
    int size = 0;
    for (it = Lin.begin(); it != Lin.end(); it++)
    {
		size += it->second.size();
    }
    for (it = Lout.begin(); it != Lout.end(); it++)
    {
		size += it->second.size();
    }
    return size;
}

/*
  Init the array stored in the DFS order, sorted by the "start" value;
  not Correct now!
*/
void Path_Hop::init_dfs_array()
{
    dfs_array = (int*) malloc(node_count * sizeof (int));
    for (int i = 0; i < node_count; i++)
    {
		dfs_array[i] = begin_2_nid(i);
    }
}

/*
  Do a batch of queries using BFS search, processing time returned;
*/
void Path_Hop::batch_query_bfs(int query_size)
{
    cout << "Generating " << query_size << " random queries...." << endl;
    srand48(time(NULL));
    int s, t;
    int left = 0;
    vector<int> src, trg;
    vector<int>::iterator sit, tit;
    int gsize = num_vertices(gra);
    while (left < query_size)
    {
		s = lrand48() % gsize;
		t = lrand48() % gsize;
		if (s == t) continue;
		src.push_back(s);
		trg.push_back(t);
		++left;
    }
    cout << "Processing queries..." << endl;
    Timer timer;
    timer.start();
    for (sit = src.begin(), tit = trg.begin(); sit != src.end(); sit++, tit++)
    {
		bfs_reach(*sit, *tit);
    }
    cout << "Query processing time (BFS) for " << query_size << " queries (ms): " << timer.stop() << endl;
}

bool Path_Hop::bfs_reach(int nid1, int nid2)
{
    if (nid1 == nid2)
		return true;
    graph_traits<Graph>::out_edge_iterator oute, oute_end;
    Vertex v = vertex(nid1, gra);
    Vertex v_next;
    std::queue<Vertex> que;
    que.push(v);
    while (!que.empty())
    {
		v = que.front();
		que.pop();
		for (tie(oute, oute_end) = out_edges(v, gra); oute != oute_end; oute++)
		{
			v_next = target(*oute, gra);
			if (gra[v_next].id == nid2)
				return true;
			que.push(v_next);
		}
    }
    return false;
}

/*
  Initialize the uncovered TC matrix, false = reachable but not yet covered, true = no needs to take care of OR reachable and is covered;
  Stop condition: all false become true;
*/
void Path_Hop::init_utc()
{
    utc = (bool**)malloc(node_count * sizeof (bool*));
    for (int i = 0; i < node_count; i++)
    {
		utc[i] = (bool*)malloc(node_count * sizeof (bool));
		for (int j = 0; j < node_count; j++)
			utc[i][j] = true;
    }
    map<int, set<int> >::iterator it_m;
    set<int>::iterator it_s;

    // construct utc according to desc-list
    for (it_m = descendants_list.begin(); it_m != descendants_list.end(); it_m++)
    {
		for (it_s = it_m->second.begin(); it_s != it_m->second.end(); it_s++)
		{
			utc[it_m->first][*it_s] = false; //false means not covered.
			tc_cnt_original++;
		}
    }
    cout << "RTC size: " << tc_cnt_original << endl;
}

/*
  Calc in_vol using the UTC matrix
*/
int Path_Hop::calc_in_vol(size_t nid)
{
    map<int, set<int> >::iterator it_m, it_m1;
    set<int>::iterator it_s, it_s1;
    it_m = ancestors_list.find(nid);
    if (it_m == ancestors_list.end())
		return 0;
    for (it_s = it_m->second.begin(); it_s != it_m->second.end(); it_s++)
    {
		if (utc[*it_s][nid])
			it_m->second.erase(it_s);
		/*else
		  {
		  it_m1 = Lout.find(*it_s);
		  if (it_m1 != Lout.end())
		  {
		  for (it_s1 = it_m1->second.begin(); it_s1 != it_m1->second.end(); it_s1++)
		  {
		  if (is_descendant(*it_s1, nid))
		  {
		  it_m->second.erase(it_s);
		  break;
		  }
		  }
		  }
		  }*/
    }
    if (it_m->second.size() <= 0)
    {
		ancestors_list.erase(it_m);
		return 0;
    }
    return it_m->second.size();
}

/*
  Calc the out_vol using the UTC matrix
*/
int Path_Hop::calc_out_vol(size_t nid)
{
    map<int, set<int> >::iterator it_m, it_m1;
    set<int>::iterator it_s, it_s1;
    it_m = descendants_list.find(nid);
    if (it_m == descendants_list.end())
		return 0;
    for (it_s = it_m->second.begin(); it_s != it_m->second.end(); it_s++)
    {
		if (utc[nid][*it_s])
			it_m->second.erase(it_s);
		/*else
		  {
		  it_m1 = Lin.find(*it_s);
		  if (it_m1 != Lin.end())
		  {
		  for (it_s1 = it_m1->second.begin(); it_s1 != it_m1->second.end(); it_s1++)
		  {
		  if (is_descendant(nid, *it_s1))
		  {
		  it_m->second.erase(it_s);
		  break;
		  }
		  }
		  }
		  }*/
    }

    if (it_m->second.size() <= 0)
    {
		descendants_list.erase(it_m);
		return 0;
    }
    return it_m->second.size();
}

/*
  Calc Density using UTC matrix;
  in_out : 1 = only in_vol, 2 = only out_vol, 3 = both in_vol and out_vol
  Here we do not count the connection from top_nid to out_vol;
*/
float Path_Hop::calc_den(size_t top_nid, size_t bottom_nid, int& in_out)
{
    map<int, set<int> >::iterator it_m1, it_m2;
    set<int>::iterator it_s1, it_s2;
    in_out = 0;
    int in_vol = calc_in_vol(top_nid);
    int out_vol = calc_out_vol(bottom_nid);
    if (in_vol > 0)
		in_out += 1;
    if (out_vol > 0)
		in_out += 2;
    if (in_out == 0)
		return 0;
    int r = 0; // the non-effective connections
    it_m1 = ancestors_list.find(top_nid);
    it_m2 = descendants_list.find(bottom_nid);
    if (it_m1 != ancestors_list.end() && it_m2 != descendants_list.end())
    {
		for (it_s1 = it_m1->second.begin(); it_s1 != it_m1->second.end(); it_s1++)
		{
			for (it_s2 = it_m2->second.begin(); it_s2 != it_m2->second.end(); it_s2++)
			{
				if (utc[*it_s1][*it_s2]) r++;
			}
		}
    }
    /*
      if (top_nid != bottom_nid && it_m2 != descendants_list.end()) // need to check only when it is a path hop
      {
      for (it_s2 = it_m2->second.begin(); it_s2 != it_m2->second.end(); it_s2++)
      {
      if (utc[top_nid][*it_s2]) r++;
      if (utc[bottom_nid][*it_s2]) r++; // already checked in calc_out_vol();
      }
      }
      int factor;
      if (top_nid == bottom_nid)
      factor = 1;
      else
      factor = 2;
    */
    //return (float) (in_vol * out_vol +in_vol + factor * out_vol - r) / (in_vol + out_vol);
    return (float) (in_vol * out_vol - r) / (in_vol + out_vol);
}

/*
  Select the densest sub bipartite graph for each hop;
  The density of such sub sgraph is returned;
*/
float Path_Hop::densest_sub_graph(int top_nid, int bottom_nid)
{
    return -1000;
    adj_nodes.clear();
    adj_nodes_revert.clear();
    adj_nodes_max.clear();
    map<int, set<int> >::iterator it_m1, it_m2, it_m;
    set<int>::iterator it_s1, it_s2, it_s;
    map<int, int>::iterator it_n;
    float den_max = -1, den = -1;
    int min_degree = INT_MAX, min_deg_node = -1, edge_cnt = 0;
    it_m1 = ancestors_list.find(top_nid);
    it_m2 = descendants_list.find(bottom_nid);
    //assert (it_m1 != ancestors_list.end());
    //assert (it_m2 != descendants_list.end());

    //initialization
    set<int> s;
    for (it_s1 = it_m1->second.begin(); it_s1 != it_m1->second.end(); it_s1++)
    {
		if (!utc[*it_s1][top_nid])
		{
			adj_nodes.insert(make_pair(*it_s1, s)); // init adjacent node list
		}
    }
    for (it_s2 = it_m2->second.begin(); it_s2 != it_m2->second.end(); it_s2++)
    {
		if (!utc[bottom_nid][*it_s2])
		{
			adj_nodes_revert.insert(make_pair(*it_s2, s)); // init adjacent node list
		}
    }
    for (it_m1 = adj_nodes.begin(); it_m1 != adj_nodes.end(); it_m1++)
    {
		for (it_m2 = adj_nodes_revert.begin(); it_m2 != adj_nodes_revert.end(); it_m2++)
		{
			if (!utc[it_m1->first][it_m2->first]) // if connection is not covered
			{
				it_m1->second.insert(it_m2->first); // add it_s2 to the adj_list of it_s1
				it_m2->second.insert(it_m1->first); // add it_s1 to the adj_list of it_s2
				edge_cnt++; // increase edge-cnt
			}
		}
    }
    while (edge_cnt > 0)
    {
		min_degree = INT_MAX;
		//find the min-degree and the first min-degree node
		for (it_m = adj_nodes.begin(); it_m != adj_nodes.end(); it_m++)
		{
			if (it_m->second.size() < min_degree)
			{
				min_degree = it_m->second.size();
				min_deg_node = it_m->first;
			}
		}
		for (it_m = adj_nodes_revert.begin(); it_m != adj_nodes_revert.end(); it_m++)
		{
			if (it_m -> second.size() < min_degree)
			{
				min_degree = it_m -> second.size();
				min_deg_node = it_m->first;
			}
		}
		//delete the min degree node
		it_m = adj_nodes.find(min_deg_node);
		if (it_m != adj_nodes.end())
		{
			edge_cnt = edge_cnt - it_m->second.size(); // update the total edge-cnt;
			for (it_s = it_m->second.begin(); it_s != it_m->second.end(); it_s++)
			{
				it_m2 = adj_nodes_revert.find(*it_s);
				it_m2->second.erase(min_deg_node);
				if (it_m2->second.size() == 0) // remove if the adj_list is empty
					adj_nodes_revert.erase(it_m2);

			}
			adj_nodes.erase(min_deg_node);
		}
		else
		{
			it_m = adj_nodes_revert.find(min_deg_node);
			edge_cnt = edge_cnt - it_m->second.size(); // update the total edge-cnt;
			for (it_s = it_m->second.begin(); it_s != it_m->second.end(); it_s++)
			{
				it_m1 = adj_nodes.find(*it_s);
				it_m1->second.erase(min_deg_node);
				if (it_m1->second.size() == 0) // remove if the adj_list is empty
					adj_nodes.erase(it_m1);
			}
			adj_nodes_revert.erase(min_deg_node);
		}
		if ((adj_nodes.size() == 0) || (adj_nodes_revert.size() == 0))
		{
			return den_max;
		}
		den = (float) edge_cnt / (adj_nodes.size() + adj_nodes_revert.size());
		if (den > den_max)
		{
			den_max = den;
			adj_nodes_max = adj_nodes;
		}
    }
    return den_max;
}

/*
 *	Fix: Merge Calc_den with Densest_sub_graph together!
 */
float Path_Hop::select_densest_hop()
{
    map<int, set<int> >::iterator it_a, it_d;
    map<int, set<nid_vol_pair, vol_Compare> >::iterator it_hd;
    set<nid_vol_pair, vol_Compare>::iterator it_s;
    int top_v, bottom_v, v;
    float den = -1, den_max = -1, den_sub = -1;
    bool is_sub_graph = false, is_sub_graph_max = false;
    int in_out = 0, in_out_max = 0; /* 1 = only in-vol, 2 = only out-vol, 3 = both in-vol and out-vol */

    //D: scan through descendants_list
    for (it_d = descendants_list.begin(); it_d != descendants_list.end(); it_d++)
    {
		//is_sub_graph = false;
		v = it_d->first;
		den = calc_den(v, v, in_out);
		/*
		if (in_out == 3)
		{
			den_sub = densest_sub_graph(v, v);
			if (den_sub > den)
			{
				den = den_sub;
				is_sub_graph = true;
			}
		}
		*/		
		if (den > den_max)
		{
			den_max = den;
			top_nid = v;
			bottom_nid = v;
			in_out_max = in_out;
			/*is_sub_graph_max = is_sub_graph;
			if (is_sub_graph_max)
			{
				adj_nodes_final = adj_nodes_max;
			}
			*/
		}
    }

    //A: scan through the ancestors_list
    for (it_a = ancestors_list.begin(); it_a != ancestors_list.end(); it_a++)
    {
		//is_sub_graph = false;
		top_v = it_a->first;
		den = calc_den(top_v, top_v, in_out);
		/*
		if (in_out == 3)
		{
			den_sub = densest_sub_graph(top_v, top_v);
			if (den_sub > den)
			{
				den = den_sub;
				is_sub_graph = true;
			}
		}
	    */   
		if (den > den_max)
		{
			den_max = den;
			top_nid = top_v;
			bottom_nid = top_v;
			in_out_max = in_out;
			/*
			is_sub_graph_max = is_sub_graph;
			if (is_sub_graph_max)
			{
				adj_nodes_final = adj_nodes_max;
			}
			*/
		}
    }
	/*
    //H: examine each hop
    //for (it_a = ancestors_list.begin(); it_a != ancestors_list.end(); it_a++)
    for (it_hd = descendants_sort_list.begin(); it_hd != descendants_sort_list.end(); it_hd++)
    {
		top_v = it_hd->first;
		//it_hd = descendants_sort_list.find(top_v);
		//if (it_hd != descendants_sort_list.end())
		//{
		//	for (it_s = it_hd->second.begin(); it_s != it_hd->second.end(); it_s++)
		//	{
		bottom_v = get_first_of_desc_sort_list(top_v);
		if (bottom_v != -1)
		{
			//is_sub_graph = false;
			den = calc_den(top_v, bottom_v, in_out);
		   
			if (in_out == 3)
			{
				den_sub = densest_sub_graph(top_v, bottom_v);
				if (den_sub > den)
				{
					den = den_sub;
					is_sub_graph = true;
				}
			}
		   
			if (den > den_max)
			{
				den_max = den;
				top_nid = top_v;
				bottom_nid = bottom_v;
				in_out_max = in_out;
				is_sub_graph_max = is_sub_graph;
				if (is_sub_graph_max)
				{
					adj_nodes_final = adj_nodes_max;
				}
			}
		}
    }
	*/
    if (in_out_max == 1)
    {
		add_to_Lout(top_nid);
		update_utc(top_nid, bottom_nid, in_out_max);
		ancestors_list.erase(top_nid);
    }
    else if (in_out_max == 2)
    {
		add_to_Lin(bottom_nid);
		update_utc(top_nid, bottom_nid, in_out_max);
		descendants_list.erase(bottom_nid);
    }
    else if (in_out_max == 3)
    {
	    
		if (is_sub_graph_max)
		{
			add_to_Lin_Lout_sub(top_nid, bottom_nid);
			update_utc_sub(top_nid, bottom_nid);
		}
		else
		{
			add_to_Lout(top_nid);
			add_to_Lin(bottom_nid);
			update_utc(top_nid, bottom_nid, in_out_max);
			ancestors_list.erase(top_nid);
			descendants_list.erase(bottom_nid);
		}
    }
	
	/*
    if (hop_cnt++ % 100 == 0)
    {
		if (top_nid == bottom_nid)
			cout << "{" << top_nid << "}";
		else
			cout << "{" << top_nid << ":" << bottom_nid << "}";
		int size = path_hop_index_size();
		cout << "Den:" << den_max << " TC:" << tc_cnt << " Sz:" << size << " R:" << (float)tc_cnt/size<< endl;
    }
	*/
    return den_max;
}

Path_Hop::Path_Hop(Graph& g) : gra(g)
{
    hop_cnt = 0;
    cout << endl << "---Path Hop Index Construction Starts---" << endl;
    cout << "Merging SCC..." << endl;
    LabelHelper helper(gra);
    helper.mergeStronglyCC();
    topo_sort();
    //reduce_graph(gra);
    cout << "Generating tree interval labeling..." << endl;
    p_interval_t = new Interval(gra);
    p_interval_nt = new Interval(gra);
    node_count = p_interval_t->createTreeLabels(false);
    p_interval_nt->createNonTreeLabels(false); //  2nd one is not correct.
    init_tree_interval();
    init_dfs_array();
    tc_cnt = tc_cnt_original = 0;
    cout << "Generating Path Hop labeling......" << endl;
    construct_descendants_list();
    construct_descendants_sort_list();
    init_utc();
}

void Path_Hop::create_label()
{
	cout << endl << "...Wait....." << endl;
    cout << "Selecting densest Path Hop...." << endl;
    while (tc_cnt < tc_cnt_original)
    {
		select_densest_hop();
    }
    remove_redundant_Lin_Lout_in_vertex();
    remove_redundant_Lin_Lout_in_tree();
    set_2_vector();
    init_parent_pointer();
    int index_size = path_hop_index_size();
    cout << endl << "---Construction finished---" << endl;
    cout << "Path Hop Label Size (no. of integers of L_in & L_out ): " << index_size << endl;
    cout << "Total Index Size ( plus 3*|V|): " << index_size + node_count * 3 << endl;
    //post_process();
}

int Path_Hop::hlss_on_utc()
{
	hl = new hlss(gra, node_count, utc);
	return hl->create_label();
}


void Path_Hop::construct_descendants_list()
{
    map<int, PairSet >::iterator it_nt;
    map<int, set<int> >::iterator it_a;
    PairSet::iterator it_s;
    set<int> descendants;
    int nid = 0;
    std::vector<Vertex>::reverse_iterator it;
    //int begin = 0, interval = 0, v = 0, node_type;
    for (it = topo.rbegin(); it != topo.rend(); ++it)
    {
		it_nt = p_interval_nt->intervals.find(gra[*it].id);
		if (it_nt->second.size() <= 0)
			continue;
		//if (has_non_tree_edge(*it)) // only for special nodes?
		{
			descendants.clear();
			for (it_s = it_nt->second.begin(); it_s != it_nt->second.end(); it_s++)
			{
				nid = interval_2_nid(it_s->first, it_s->second);
				descendants.insert(nid);
				//add to ancestors_list
				it_a = ancestors_list.find(nid);
				if (it_a == ancestors_list.end())
				{
					set<int> s;
					s.insert(it_nt->first);
					ancestors_list.insert(make_pair(nid, s));
				}
				else
				{
					it_a->second.insert(it_nt->first);
				}
			}
			if (descendants.size() > 0)
				descendants_list.insert(make_pair(it_nt->first, descendants));
		}
    }
}

void Path_Hop::update_utc_sub(int top_nid, int bottom_nid)
{
    //cout << "update utc sub:" << top_nid << ":" << bottom_nid << endl;
    map<int, set<int> >::iterator it_m, it_m1;
    set<int>::iterator it_s;
    int x, y;
    for (it_m = adj_nodes_final.begin(); it_m != adj_nodes_final.end(); it_m++)
    {
		x = it_m->first;
		if (!utc[x][top_nid])
		{
			utc[x][top_nid] = true;
			tc_cnt++;
		}

		for (it_s = it_m->second.begin(); it_s != it_m->second.end(); it_s++)
		{
			y = *it_s;
			if (!utc[x][y])
			{
				utc[x][y] = true;
				tc_cnt++;
			}
			if (!utc[top_nid][y])
			{
				utc[top_nid][y] = true;
				tc_cnt++;
			}
			if (!utc[bottom_nid][y])
			{
				utc[bottom_nid][y] = true;
				tc_cnt++;
			}
		}
    }
}


/*
  Update the utc matrix, convert these connections from false to true (processed):
  1. in_vol  --> top_nid;
  2. bottom_nid --> out_vol;
  3. in_vol  --> out_vol;
  4. top_nid --> out_vol;
*/
void Path_Hop::update_utc(int top_nid, int bottom_nid, int hop_type)
{

    map<int, set<int> >::iterator it_m1, it_m2;
    set<int>::iterator it_s1, it_s2;
    it_m1 = ancestors_list.find(top_nid);
    it_m2 = descendants_list.find(bottom_nid);
    if (hop_type == 1)
    {
		if (it_m1 != ancestors_list.end())
		{
			for (it_s1 = it_m1->second.begin(); it_s1 != it_m1->second.end(); it_s1++)
			{
				if (!utc[*it_s1][top_nid])
				{
					utc[*it_s1][top_nid] = true;
					tc_cnt++;
				}
			}
		}
    }
    else if (hop_type == 2)
    {
		if (it_m2 != descendants_list.end())
		{
			for (it_s2 = it_m2->second.begin(); it_s2 != it_m2->second.end(); it_s2++)
			{
				if (!utc[bottom_nid][*it_s2])
				{
					utc[bottom_nid][*it_s2] = true;
					tc_cnt++;
				}
			}
		}
    }
    else if (hop_type == 3)
    {
		//if (it_m1 == ancestors_list.end()) return;
		//if (it_m2 == descendants_list.end()) return;
		if ((it_m1 != ancestors_list.end()) && (it_m2 != descendants_list.end()))
		{
			for (it_s1 = it_m1->second.begin(); it_s1 != it_m1->second.end(); it_s1++)
			{
				if (!utc[*it_s1][top_nid])
				{
					utc[*it_s1][top_nid] = true;
					tc_cnt++;
				}
			}
			for (it_s2 = it_m2->second.begin(); it_s2 != it_m2->second.end(); it_s2++)
			{
				if (!utc[bottom_nid][*it_s2])
				{
					utc[bottom_nid][*it_s2] = true;
					tc_cnt++;
				}
				if (!utc[top_nid][*it_s2])
				{
					utc[top_nid][*it_s2] = true;
					tc_cnt++;
				}
			}
			for (it_s1 = it_m1->second.begin(); it_s1 != it_m1->second.end(); it_s1++)
			{
				for (it_s2 = it_m2->second.begin(); it_s2 != it_m2->second.end(); it_s2++)
				{
					if (!utc[*it_s1][*it_s2])
					{
						utc[*it_s1][*it_s2] = true;
						tc_cnt++;
					}
				}
			}
		}
    }
}

void Path_Hop::add_to_Lin_Lout_sub(int top_nid, int bottom_nid)
{
    map<int, set<int> >::iterator it_m, it;
    set<int>::iterator it_s, it_s2;
    set<int> X, Y;
    for (it_m = adj_nodes_final.begin(); it_m != adj_nodes_final.end(); it_m++)
    {
		X.insert(it_m->first);
		for (it_s = it_m->second.begin(); it_s != it_m->second.end(); it_s++)
		{
			Y.insert(*it_s);
		}
    }

    for (it_s = X.begin(); it_s != X.end(); it_s++)
    {
		if (utc[*it_s][top_nid]) continue;
		it = Lout.find(*it_s);
		if (it == Lout.end())
		{
			set<int> s;
			s.insert(top_nid);
			Lout.insert(make_pair(*it_s, s));
		}
		else
		{
			//check if *it_m could reach any anscestors of v by non-tree egdes; If yes, then no need to add v any more;
			bool already_reachable = false;
			for (it_s2 = it->second.begin(); it_s2 != it->second.end(); it_s2++)
			{
				if (is_descendant(*it_s2, top_nid))
				{
					already_reachable = true;
					break;
				}
			}
			if (!already_reachable)
			{
				it->second.insert(top_nid);
			}
		}
    }
    for (it_s = Y.begin(); it_s != Y.end(); it_s++)
    {
		if (utc[bottom_nid][*it_s]) continue;
		it = Lin.find(*it_s);
		if (it == Lin.end())
		{
			set<int> s;
			s.insert(bottom_nid);
			Lin.insert(make_pair(*it_s, s));
		}
		else
		{
			//check if *it_m could reach any anscestors of v by non-tree egdes; If yes, then no need to add v any more;
			bool already_reachable = false;
			for (it_s2 = it->second.begin(); it_s2 != it->second.end(); it_s2++)
			{
				if (is_descendant(bottom_nid, *it_s2))
				{
					already_reachable = true;
					break;
				}
			}
			if (!already_reachable)
			{
				it->second.insert(bottom_nid);
			}
		}
    }
}


/*
  add v to Lin of corresponding nodes
*/
void Path_Hop::add_to_Lin(int v)
{
    int begin = 0, interval = 0, n = 0, node_type = 0;
    map<int, set<int> >::iterator it_m, it;
    set<int>::iterator it_s, it_s2;
    it_m = descendants_list.find(v);
    if (it_m == descendants_list.end())
		return;
    bool already_reachable;
    for (it_s = it_m->second.begin(); it_s != it_m->second.end(); it_s++)
    {
		if (utc[v][*it_s]) continue;
		it = Lin.find(*it_s);
		if (it == Lin.end())
		{
			set<int> s;
			s.insert(v);
			Lin.insert(make_pair(*it_s, s));
		}
		else
		{
			already_reachable = false;
			for (it_s2 = it->second.begin(); it_s2 != it->second.end(); it_s2++)
			{
				//assert(v != *it_s2);
				if (is_descendant(v, *it_s2))
				{
					already_reachable = true;
					break;
				}
			}
			if (!already_reachable)
			{
				it->second.insert(v);
			}
		}

    }
}

/* 
   add v to Lout of corresponding nodes
*/
void Path_Hop::add_to_Lout(int v)
{
    map<int, set<int> >::iterator it_m, it;
    set<int>::iterator it_s, it_s2;
    it_m = ancestors_list.find(v);
    if (it_m == ancestors_list.end())
		return;
    bool already_reachable;
    for (it_s = it_m->second.begin(); it_s != it_m->second.end(); it_s++)
    {
		if (utc[*it_s][v]) continue;
		it = Lout.find(*it_s);
		if (it == Lout.end())
		{
			set<int> s;
			s.insert(v);
			Lout.insert(make_pair(*it_s, s));
		}
		else
		{
			//currently we do not do this optimization
			//check if *it_m could reach any anscestors of v by non-tree egdes; If yes, then no need to add v any more;
			already_reachable = false;
			for (it_s2 = it->second.begin(); it_s2 != it->second.end(); it_s2++)
			{
				if (is_descendant(*it_s2, v))
				{
					already_reachable = true;
					break;
				}
			}
			if (!already_reachable)
			{
				it->second.insert(v);
			}
		}
    }
}

void Path_Hop::user_query()
{
    while (true)
    {
		int nid1, nid2;
		cout << "Pls input two nodes' id (enter two -1 to stop):" << endl;
		cin >> nid1 >> nid2;
		if (nid1 == -1) return;
		cout << "Your input is: " << nid1 << "," << nid2;
		if (reach(nid1, nid2))
			cout << " ---> Yes!" << endl;
		else
			cout << " ---> No!" << endl;
    }
}

int Path_Hop::check_correctness(int test_size)
{
    vector<Vertex> nodes;
    int gsize = num_vertices(gra);
    // get the graph size
    graph_traits<Graph>::vertex_iterator vi, vi_end;
    for (tie(vi, vi_end) = vertices(gra); vi != vi_end; ++vi)
    {
		nodes.push_back(*vi);
		// store each vertex into a vector
    }
    int s, t;
    // source and targetindex
    int left = test_size;
    int err = 0;
    //Interval intl(gra);
    //intl.createLabels(false);
    while (left > 0)
    {
		// randomly genera a query
		s = rand() % gsize;
		t = rand() % gsize;
		if (s == t)
		{
			// not to test the same node
			continue;
		}
		if (p_intl->reachable(nodes[s], nodes[t]) != this->reach(gra[nodes[s]].id, gra[nodes[t]].id))
		{
			err++;
			cout << err << ":  " << gra[nodes[s]].id << " (" << gra[nodes[s]].begin << "," << gra[nodes[s]].end << ") "
				 << "--->" << gra[nodes[t]].id << " (" << gra[nodes[t]].begin << "," << gra[nodes[t]].end << ") "
				 << "  [" << p_intl->reachable(nodes[s], nodes[t]) << ", " << reach(gra[nodes[s]].id, gra[nodes[t]].id) << "]" << endl;
		}
		left--;
    }
    return err;
}

/*
  Do a batch of queries upon path-hop labels, processing time returned;
*/
void Path_Hop::batch_query(int query_size)
{
    cout << "Generating " << query_size << " random queries...." << endl;
    srand48(time(NULL));
    int s, t;
    int left = 0;
    vector<int> src, trg;
    vector<int>::iterator sit, tit;
    int gsize = num_vertices(gra);
    while (left < query_size)
    {
		s = lrand48() % gsize;
		t = lrand48() % gsize;
		if (s == t) continue;
		src.push_back(s);
		trg.push_back(t);
		++left;
    }
    cout << "Processing queries..." << endl;
    parent_jump = 0;
    Timer timer;
	timer.start();
    for (sit = src.begin(), tit = trg.begin(); sit != src.end(); sit++, tit++)
    {
		reach(*sit, *tit);
    }
	cout << endl << "Query processing time: " << timer.stop() << endl;
}

/*
  Test if nid1 reaches nid2, Problem fix: too slow!
*/
bool Path_Hop::reach(int nid1, int nid2)
{
    if (tree_interval[nid1][0] <= tree_interval[nid2][0] && tree_interval[nid1][1] >= tree_interval[nid2][1]) return true;
    vector<int>::iterator p1, p2;
    ancestors_Lin_vec.clear();
    descendants_Lout_vec.clear();
    descendants_Lout_vec.push_back(nid1);
    descendants_Lout_vec.insert(descendants_Lout_vec.begin(), Lout_vec[nid1].begin(), Lout_vec[nid1].end());
    ancestors_Lin_vec.push_back(nid2);
    ancestors_Lin_vec.insert(ancestors_Lin_vec.begin(), Lin_vec[nid2].begin(), Lin_vec[nid2].end());
    int parent = tree_interval[nid2][2];
    while (parent != -1)
    {
		ancestors_Lin_vec.insert(ancestors_Lin_vec.begin(), Lin_vec[parent].begin(), Lin_vec[parent].end());
		parent = tree_interval[parent][2];
		parent_jump++;
    }
    int start, end;
    for (p1 = descendants_Lout_vec.begin(); p1 != descendants_Lout_vec.end(); p1++)
    {
		start = tree_interval[*p1][0];
		end = tree_interval[*p1][1];
		for (p2 = ancestors_Lin_vec.begin(); p2 != ancestors_Lin_vec.end(); p2++)
		{
			if (start <= tree_interval[*p2][0] && end >= tree_interval[*p2][1]) return true;
		}
    }
    return false;
}
